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We consider solitary wave solutions to the Dirac-Coulomb system both from physical and mathematical points 
of view. Fermions interacting with gravity in the Newtonian limit are described by the model of Dirac fermions 
with the Coulomb attraction. This model also appears in certain condensed matter systems with emergent Dirac 
fermions interacting via optical phonons. In this model, the classical soliton solutions of equations of motion 
describe the physical objects that may be called polarons in analogy to the solutions of the Choquard-Pekar 
equation. We develop analytical methods for the Dirac-Coulomb system, showing that the no-node gap solitons 
for sufficiently small values of charge are linearly (or spectrally) stable. 

PACS numbers: 



I. INTRODUCTION 

The Dirac -Maxwell system, as well as models of self-interacting fermion fields such as the massive Thirring model 
ll49ll and the Soler model j4(J, have been attracting the interest of both physicists and mathematicians for many years. 
In these models, just like in the nonlinear Schrodinger equation, there exist localized solitary wave solutions of the 
form (f)(x)e~ tut , with (f>(x) exponentially localized in space. These solutions are due to the nonlinear nature of the 
systems and may encode certain properties of the theory which can not be obtained via perturbative analysis. The role 
of these classical solutions in the High Energy Physics, as models of extended fermions and confined Dirac quarks, 
has long been the topic of intense discussion, going back to iB 12811321 l40l l4ll l45ll . 

In the nonlinear Dirac equation (Soler model I46I0 . 

it = -ia ■ VC + (m - CO/30 

with C(x,t) G C 4 and x G M 3 , the existence of solitary wave solutions was shown in ITii 14^1 . The question of 
stability of these solitary waves attracted much attention of physicists for many years, but only partial results were 
obtained. (Let us mention that, in contrast, the questions on linear, orbital, and even asymptotic stability of solitary 
waves in the nonlinear Schrodinger equation are essentially settled lfl2l I2TI l27i l50ll .) The minimum e nerg y-charge 
solitary wave corresponding to w c = 0.936m was originally proposed as a model of extended fermions fl46|], but was 
later shown to be unstable @|. (General results on instability of minimal energy solitary waves are in jlStl .) The 
analysis of stability of Dirac solitary waves with respect to particular families of perturbations (such as dilations) has 
been performed in HI H3l ; yet, neither the linear stability (spectrum of the equation linearized at a soliton) nor 
orbital nor asymptotic stability have been understood. Let us also mention the related numerical results 1I3TI5I [l6ll . 

Using the approach described in the second part of the present paper (Section IV BK one expects to prove linear 
stability of small solitary wave solutions to the nonlinear Dirac equation in ID (the massive Gross-Neveu model I32T0 . 
For the nonlinear Dirac equation in 3D, the small charge solitary waves (w < m) are linearly unstable (this can be 
shown using the methods of lfl8ll ). but it is reasonable to expect that the solitary waves with oj < uj c = 0.936m are 
linearly stable: as u> decreases below w c , one positive and one negative eigenvalue of the linearized operator collide at 
A = and then re-emerge as two purely imaginary eigenvalues lfl7l fl9"tl . 

Let us mention related results on the Einstein-Dirac equations. This system was shown to have particle-like solutions 
which are linearly stable with respect to spherically symmetric perturbations ll26ll . 

For the Dirac-Maxwell system, the localized solutions with ui G (— m, m) have been shown to exist, first numer- 
ically and then analytically, in 01 [24], HH HH] . Some of the solitary waves in the Dirac-Maxwell system seem to be 
linearly stable, at least for small values of the charge, when u > — m; the approach developed in Section IVBl is ap- 
plicable to this system, yielding this result. Yet, it seems that there is no direct physical meaning of the solitary waves 
in the classical Dirac-Maxwell system. These localized states are formed due to attraction of the spinor field to itself, 
which takes place at the energies ui > —m. At the same time, at energies ui < in the quantum theory the antiparti- 
cles would be appearing. To describe antiparticles, one needs to change the order of the fermion creation-annihilation 
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operators. This results in the additional change of sign at the scalar potential. Thus, in the quantum theory, instead of 
self-attraction which could lead to the creation of a certain localized mode, one again ends up with the self -repulsion 
in the antiparticle sector. This anticommuting nature of fermion variables is completely disregarded in the classical 
Dirac-Maxwell system. 

Let us also mention the positronium, which is the bound state of electron and positron in the quantum Dirac-Maxwell 
system. This state has nothing to do with the solitary waves discussed in this paper. Positronium can be described 
perturbatively, by the Dirac equation in external Coulomb field; this is in contrast to the solitary wave solutions, which 
are described by Dirac equation in the potential created by the spinor field itself. The nature of the solitary waves is, 
therefore, essentially nonlinear, unlike the eigenstate of the Dirac equation corresponding to positronium. 

In the present paper, we make an attempt to find out the place in the quantum field theory where the classical- 
theoretical localized solutions could play a role in the dynamics. We focus on the fermion system in which massive 
electrons attract each other via the Coulomb-like interaction. This toy model may play a role in several physical 
problems to be discussed briefly below. 

It is well known that the relativistic Dirac fermions may emerge in the condensed matter systems. This occurs, 
for example, on the boundary of the 3D topological insulators and in some two-dimensional structures like graphene 
lEM l36ll . Moreover, massless Dirac fermions appear in 3D materials (see, for example, |p4fl and references therein) in 
the intermediate state at the phase transition between a topological and a normal insulator. Recently, the existence of 
massless fermions at the phase transition between the insulator states with different values of the topological invariants 
has been proved within a wide class of relativistic models 115 511 . When the interaction with other fields is taken into 
account, these massless particles gain the mass. Fermion excitations in various 3D materials and other systems interact 
with phonons l3lll ; as a result, the attractive interaction between the fermions appear. Such an interaction gives rise to 
the formation of the Cooper pairs in the microscopic mechanism of the superconductivity ll3lll . The interaction may 
also have a form similar to the Coulomb attraction (for example, in the case of the exchange by optical phonons lfl5lo . 
This interaction for the case of the conventional nonrelativistic electron leads, in particular, to the formation of the 
polaron JUS- 

In this paper, we show that the interaction of fermions with optical phonons (polaron problem for Dirac field) in 
a certain approximation is equivalent to the classical equations of motion for the massive Dirac spinor field with 
Coulomb attraction (see also description of Landau-Pekar approach in l30ll38ll ). 

In the above model, the Lorentz symmetry is broken due to the Coulomb forces. The relativistically-invariant model 
of a similar type is given by the Dirac fermions interacting with the gravitational field. This problem may be thought 
of as a true relativistic polaron problem. We show that in the Newtonian limit we again arrive at the system of Dirac 
fermions interacting via Coulomb-like forces. Such polarons may play the role of nonperturbative solutions in certain 
unified theories. 

We show that the Dirac equation with the Coulomb attraction may emerge in the semiclassical description of the 
corresponding quantum problems. 

In the second part of the paper, we analyze the stability of solitary wave solutions in the resulting Dirac-Coulomb 
system, and give the analytic justification of their linear stability. Namely, we show that certain solitary waves are 
linearly stable (or spectrally stable): the spectrum of the equation linearized at a particular solitary wave has no 
eigenvalues with positive real part. Our approach is based on the fact that the nonrelativistic limit of the Dirac-Coulomb 
system is the Choquard-Pekar equation, and in particular the solitary wave solutions to the Dirac-Coulomb are obtained 
as a bifurcation from the solitary waves of the Choquard-Pekar equation. (It is worth mentioning that the Choquard- 
Pekar equation appears in the conventional polaron problem ll30ll3~8ll .) It also follows that in the nonrelativistic limit the 
families of eigenvalues of Dirac-Coulomb linearized at a soliton are deformations of eigenvalue families corresponding 
to Choquard-Pekar. The latter could be analyzed via the Vakhitov-Kolokolov stability criterion lf30tl - The delicate part 
is the absence of bifurcations of eigenvalues from the continuous spectrum; we explain that this follows from the 
limiting absorption principle for the free Dirac operator. 

We emphasize that this is the first result for the linear stability of spatially localized fermion modes. We expect 
that this approach is applicable to other systems, such as the nonlinear Dirac equation and the Dirac-Maxwell system. 
Let us also mention recent results on asymptotic stability of solitary wave solutions to the nonlinear Dirac equation in 
ID and in 3D lflol[39ll . proved under the assumption that the corresponding solitary waves are linearly stable. Other 
related results on linear stability and instability for the nonlinear Dirac equation are in lp7L Fnl [Til I20I1 . 

The paper is organized as follows. In Section HH we briefly describe the model under investigation. In Section Hill 
we discuss the appearance of the considered solitary waves in the problem of quantum Dirac fermion interacting with 
optical phonons. In Section[lV]we consider the relation of solitary waves to the gravitational polaron. In SectionlVlwe 
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address the question of stability of these solitary waves. The conclusions are in Section[VT] In Appendix, we give the 
details of the Vakhitov-Kolokolov stability criterion f50ll in its application to the Choquard-Pekar equation. 



II. DIRAC-COULOMB SYSTEM 

Below, we choose the units so that K = c = 1. We consider the model with the action 



S E = fd 3 xdtC(iD v -m)C-^ fd 3 x dt ^T(d 3 Lp) 2 , (II. 1) 

J J J= i 



where Dirac fermion field interacts with itself via an instantaneous Coulomb interaction, x G R 3 , t G R, e 2 is the 
coupling constant, ((x, t) G C 4 is a four-component Dirac field, ip(x, t) G R is the Coulomb field, and 

3 
3=1 

where da — Jj, dj — -^j, 1 < j < 3; the Dirac matrices 7^, < fx < 3, satisfy the Euclidean-Clifford algebra 
{7^1 7^} = 2g M,y , with g MI/ the inverse of the metric tensor g^,, = diag[l, —1, —1, —1]. It is worth mentioning that in 
this model the Lorentz symmetry is broken due to the Coulomb forces. 

The dynamical equations corresponding to (III, lb are given by the following Dirac-Coulomb system: 



idtC = -ia ■ VC + mpQ + etp(, 

= e(*C, ( ' 

where ((x,t) G C 4 , ip(x,t) G R, x G R 3 . Above, e is the charge of the spinor field, <p(x,t) is the (real-valued) 
external scalar field (such as the potential of the electric field in R 3 ). The self-adjoint Dirac matrices aP and (3 are 
related to 7 M by 

7 ^= 7 V, 1 < j < 3; 7°=^. 

They satisfy (a J ) 2 = (3 2 = 7 4 , a> a k + a k a^ = 2I 4 S jk , a 3 (5 + Pa j = 0; 1 < j, k < 3. Above, C = (/?C)* = CP, 
with C* the hermitian conjugate of £. A particular choice of the Dirac matrices does not matter; we take the Dirac 
matrices in the common form 



crA h 

*i 0)> P {0 -h 

r ■ n n u v /0 l\ A) -A A 

where J2 is a 2 x 2 unit matrix and er, are the Pauli matrices: o"i = „ , 01 = . _ , (T3 



(II.3) 



IIL POLARONS FORMED DUE TO INTERACTION OF FERMION FIELD WITH OPTICAL PHONONS 

A. Field-theoretical description of the generalized Frohlich model 

It was mentioned in the introduction that Dirac fermions may emerge in various 3D systems at the phase transition 
between the insulating states with different values of momentum space topological invariants. Similarly to the ordinary 
electrons in crystals, these Dirac fermions may interact with optical phonons thus giving rise to the polaron problem. 
The model Hamiltonian for a Dirac particle interacting with optical phonons can be obtained via the generalization of 
the conventional Frohlich hamiltonian 12311: 



W = EEs + [7VPi+m7%+ (HIT) 
i=i v 



k,p 1 
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Here a and il are coupling constants, cj. are the electron creation operators, and a k are the phonon creation operators. 
We introduce the phonon field tp(x) = iV2D,^2 k j^a,ke lkx + H.C. and the electron field ip{x) = ^2 p c p e lkx . 

We express Tre~ lWT , where T is time, as a functional integral. This is done as follows. First, we subdivide 
the time interval [0, T] into smaller intervals AT and represent e~ inT = e ~' MAT . . . e _lWAT . Next, we substitute 
1 = fdridf]d(d^e~''^~ f ' r] \r]0(r](\, where the coherent states are \r)Q — e vc +< = a |0), ( G C, and r] is the 
Grassmann variable. Functional integral over £, r] appears. Then, the integration variables ip and ip are introduced as a 
result of the action of the corresponding operator fields on the coherent states. Therefore, Frohlich hamiltonian gives 
rise to the quantum field theory of interacting fermion and phonon fields, with the partition function given by 



Z = Jdipd-tpdipexp(i Jd 3 xdtiP[i$-m-e- i p]i>+ % - Jd 3 xdt[-^(\7p>) 2 - (V^) 2 ]), (III.2) 

where = 7^9,,, with the summation over \i = 0, . . . , 3; (V^) 2 = X^=i( Jf - ) 2 - Above, we denote e = a^/^j- 
We formalize this as follows: 

Lemma III.l. The quantum-mechanical system with the Frohlich hamiltonian iIH.lt is equivalent to the field theory 
with the partition function HIH.2\) . 

It is worth mentioning that usually the nontrivial vierbein appears when Dirac fermions emerge in condensed matter 
models ll5lll . Moreover, this vierbein fluctuates. Under certain circumstances, such fluctuations may give rise to 
emergent gravity ll5lll . We, therefore, imply that the emergent vierbein has small fluctuations and can be transferred to 
unity matrix via the rescaling of space and time coordinates. In this model the emergent Lorentz symmetry is broken 
due to the Coulomb forces. That is why we deal with an exotic situation: the fermion hamiltonian is of the Dirac form, 
while the interactions are purely non-relativistic. 

In the low energy approximation E <C fi, we arrive at the partition function with the action given by Eq. dll.lt : 



Z = j dtpdipdtpexp(i jd 3 xdt^[i0 -m- e 7 V]^ - | j d 3 x dt [(V ip) 2 }^. (111.3) 
After the Wick rotation (t — > —it, ip — > ip) to Euclidean space-time we arrive at 



Z = 



d^) dtp d^exp( - Jd 3 xdt ip([T°(d + iep) + T j d ] + m)V> + ^ Jd 3 x dt {Vpf^j . {Ill A) 



Here the Euclidean gamma-matrices T M satisfy {r M , T"} = 26^", < n, v < 3. 

Remark III. 2. In fact, the integration in dHI.41 > is not convergent due to the positive sign at the kinetic term for p. 
This is exactly the same problem as for the Euclidean functional integral for the gravitational theory with Einstein- 
Hilbert action (see below). This shows that the theory defined by the partition function JIII.3I > can be considered 
only as an effective low energy model. At some (large) energies, the action has to be redefined in order to make the 
Euclidean functional integral convergent. As well as for the quantum gravity, this can be done if the term with higher 
powers of dp is added to the action. This, in turn, regularizes the Coulomb interaction at small distances. In physical 
applications, such additional term in the action appears at scales at which the attraction due to optical phonons no 
longer dominates, and some other interactions come into play (say, the Coulomb repulsion due to photons). 

B. Application of semiclassical methods to the model 

Here it is important that we consider the phonon field constant in time. We come to the model with the following 
partition function: 

Z = Jd^d^df exp (i J2 V T Jd 3 x V>+ [j? - H^ipr, 

-yd 3 xdt[(\'p) 2 }), (111.5) 
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where 

= 7°(7jO | i + m + e 7 V)- 

Here the system is considered with boundary conditions anti-periodic in time: ip(t + T,x) = —ip(t, x). We use the 
decomposition tp(t 7 x) = Y^ v =§ (2fc+i) e^ tnt ip v (x). We represent -0 as ^0*0 = ^„ °ri,n$n{ x )> where is the 
eigenfunction of "H^, corresponding to the eigenvalue E£ and normalized to unity ( Jd 3 x ^^^n = 1): 

f _ » E Tc,,»[>)--B*]c,,„-| /A * [(V V ) 2 ] 

Z = dcdcdipe v >" 



Integrating out Grassmann variables c n we come to: 

Z = Jdipe-if d3xdt ^ v ^ 2 %, n [r]-E^} (III.6) 



const /^e-i/ d3 * dt K v ^n„e-^ T co S [ T 



2 J ' 

where ilo depends on the details of the regularization but does not depend neither on T nor on the spectrum in 
continuum limit. The values £ n depend on the parameters of the hamiltonian, with the index n enumerating these 
values. 

Eq. dHI.6b is proved as follows. Remind that r\ = £ (2k + 1). The product over k can be calculated as follows: 



Hk=-N, jv[^(2fc + 1) + = n^.jy, N -(2k + i) n fc [i 



T v , - nJ v ,v T v /V 7r(2fc + l) 



v n^L_ JV+1 (2fc + 1) j cos(^T/2) a (-j e"^ v cos(^T/2), (III.7) 

where T = 2Na, and a is the lattice spacing. Thus, 

Det(i<9 - Hp) = e- n ° T Il n cos(E%T/2), (III.8) 

We get (see also JH El): 

Z = const e-° oT ^ fdcpexpf— fd 3 x[(Wip) 2 ] + ^J2 E n- iT zZ K ^ E n)- ( m ' 9 ) 

{if„}=0,l n n 

Following JH, we interpret Eq. ( |III.9t as follows. K n represents the number of occupied states with the energy 
E%. These numbers may be or 1. The term J2 n E% vanishes if ip = 0, since in this case the values E n come in pairs 
with opposite signs. 

Suppose that we can evaluate integral over <p in a semiclassical way. This can be done, for example, in the weak 

2 

coupling approximation when |- <C 1. Then we come to the following variational problem: 

= S{ - \d*x + £ (f - KJ%) } = S fd 3 x{ -Zf-Z (K n £U v C n + 

In the right-hand side, the variation over £„ with the constraint Jd 3 x CnCn — 1 gives* m addition, the one-fermion 
wave functions £„. The variational problem can be written as 



= 8 Jd 3 x { £(i - K n )(+{H V - A„]Cn - ^y^} = 8 Jd 3 x { ]T(A„ - 1)^0 -m- e 7 V]C n 
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Here X n are the Lagrange multipliers. We introduced the dependence on time into £: the variation is performed with 
respect to the functions of the form f„ = e~ lXnt ((x) and with respect to the time-independent phonon cloud ip. In 
this form, the functional to be used in the variational problem almost coincides with the action from (IIII.2I ). The 
difference is that we assume the special form of £ ~ e~ lXt and also that ip does not depend on time. Also, instead of 
the Grassmann variables, we substitute ordinary wave functions and take into account filling factors for the fermion 
states. 

The additional constraint is that the wave functions („ are different, so that there are no states that are occupied more 
than once. The variation is performed with the fixed numbers K n . After the variational problem is solved one says 
that the state with the found wave function £„ is occupied if K n 7^ for the corresponding value of n. In Eq. ( 1HI.9I I 
we need to sum up all such configurations with different arrays K n . The calculated values of ip are to be substituted 
into the exponent in Eq. dHI.91 l while the values of E n are given by E n = This is the kind of the Hartree - 

Fock approximation for the given model. 

We come to the following result: 

Theorem III.l. The partition function for the system of Dirac electrons interacting with optical phonons at low 
energies E <C £1 is given by Eq. MIL 9h In the weak coupling, the integral over ip in ( 1/77.91 ) is taken in the stationary 
phase approximation that leads to the variational problem Eq. MII.KM . 

Remark III. 3. CP-invariance implies that for any state n with the energy E% there exists a state h with the energy 
E'f = —E~ v . In particular, for (p = we obtain J2 n = 0. States with positive E n are interpreted as electrons, 
while the states with negative E n correspond to holes. The vacuum state is in this case the state with all negative 
levels E n occupied. The situation is changed when <p ^ 0. However, if maxe<^ < m, the states with E„ > are also 
interpreted as electrons while the states with E n < are interpreted as holes. If max e<p > 2m, then there could be 
states which can not be considered as either electrons or holes. Instead, these states correspond to the Schwinger pair 
creation process. The appearance of such states, however, may be avoided in the weak coupling, when e is small and, 
therefore, almost always eip < m. That is why in the weak coupling the vacuum can again be considered as the state 

with all negative levels of E n occupied and all positive levels of E n empty. At large enough values of a = this 
pattern may be changed due to the creation of pairs that may lead to the change of vacuum. In this case the fermion 
condensate may appear; the description of the theory in terms of the collection of one-fermion states is no longer 
relevant. 



Now let us consider the usual polaron problem, i.e. the problem of one electron interacting with the phonon cloud. 
Recall that Z = Tre"^. Therefore, sum in Eq. flIII.9t corresponds to the sum over many-fermion states. The state 
with A',j ac = i(l — s ig' n -En) corresponds to the vacuum. Then the state with K n {q) = K^ ac + 8 nq for some E q > 
corresponds to the state that consists of vacuum (Dirac sea of negative levels) and the bound state of one electron and 
the phonon cloud surrounding it. This is a polaron. Vacuum is translational and CP-invariant. That is why <p vac = 0. 
The vacuum energy is defined as E vac = ^2 n K„ ac E®. The polaron energy is equal to 



where ip q is defined by solving variational problem Eq. ( IIII- 1 Ob - Infinite vacuum energy has to be subtracted from E q . 
The quantity £ q = E q — E vac is finite and is considered as a renormalized polaron state energy in physical applications. 
In the general case, the renormalized polaron energy contains the contribution from the virtual electron-hole pairs that 
are born when large enough ip q appears. Also, the interaction with the Dirac sea gives contribution to £ q . However, 
under certain circumstances, these contributions can be neglected. Roughly speaking, this happens in the so-called 
quenched approximation when the interactions between different fermions are neglected while the interaction between 
phonon cloud and the electron is taken into account. In this case, the polaron energy is calculated as 



C. One-polaron problem 




(III. 10) 




(III. 11) 



where ip q is calculated via the variational problem 




(III. 12) 
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This problem, in turn, leads to the following equations: 

idtC = -iot ■ VC + mfiC,+ etp q (x,t)((x,i), (111.13) 
where ifl 3 , t) £ C 4 . The potential ip q (x, t) g R is generated by the spinor field itself: 

Aip q (x,t) = e|C(M)| 2 , seM 3 , (111.14) 

where A = X) J= i ^ s me Laplace operator in M 3 . 

That is why we come to Eq. ( III. 2b with the important constraint on the wave function (: 

Jd 3 xQ + C, = l. (111.15) 

We arrive at the following statement: 

Lemma III.4. Let us consider the problem of bound states of one electron surrounded by the phonon cloud (polaron) 
in the field theory with partition function MII.2h In the low energy approximation E <C f2 the system of equations 
MIL 1 51 ), UH.14t , MIL 751 ) as weW as variational problem MII.12i . solve this problem in the first order approximation 
of the weak coupling expansion. 

Remark 111.5. It is important that only those solutions of the original system of equations Eq. ( III. 2b that are normalized 
according to Eq. ( IIII. 1 5b have the physical meaning. 

It is worth mentioning that the given variational problem may be relevant for the solution of polaron problem not 
only in the weak coupling regime (see, for exampl e. ll23L f30l L38h '). In non-relativistic case this variational problem has 
appeared in the approach due to Landau and Pekar I30ll38ll . However, in the original papers lf30l L38h the trial functions 
were used for the minimization, while we consider the given variational problem exactly. 

It is instructive to consider how Eq. ( IIII. 12b appears from the consideration of the two-point Green function 



G(t a -ti) = i Jdi,d^d^e^ d3xdt y' [l0 - m -^°^- i ^}^ + {t 1 ,x)i:{t 2 ,x)d 3 x. (01.16) 

It is worth mentioning that the consideration of arbitrary values of t\ , t% requires more complicated technique, when tp 
is constant as a function of time except for the points t\, t-x. This is because the insertion of ip and xp disturbs vacuum 
in such a way that the value of ip is changed at t\ and t 2 . This technique uses the so-called Floquet indices and will 
be applied in the next section to the consideration of gravitational polarons, which are the relativistic generalizations 
of the objects considered in this section. Here we restrict ourselves by the consideration of the case t\ = 0, ti = T. 
Then 

OCT) = C ° DS l e r T E /#e--/^^e-^(--»)^ £ Z£_!J. (HU7) 
il Z L — ' / 77 — Ljq 

{A'„}=0,1 7j=f (2fc+l),g ' q 



Using Poisson summation formula we get 



r;=f (2fe+l) 



-iTe-™^ 



V-E$ 1 



Here it is implied that energy levels have small imaginary parts (as usual in quantum field theory). After the Wick 
rotation (T — > —i/T,T is temperature) the given expression would become the usual finite temperature Matsubara 
Green function. Therefore, we arrive at 



Y. fdf e lT (-^ x k«K-e$) 

G(T) = ^ E Sd(p ^T(-f^ X ^. + ^ n Sf-^ n K n ES) 

{K n }=0,l 



(III. 18) 
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Here in each term of the summation over K n the field ip is to be determined in the stationary phase approximation via 
solving the variational problem dill- 1 Ob . Again, in the quenched approximation we come to 

G(T) = W tbpe-^V** = J2 Z q e- i£ < T . 

q q 

Here £ q is the energy of the polaron in the q-th state calculated according to Lemma IIII.4I Factors Z q are the pre- 
exponential factors of the stationary phase approximation. 



IV. GRAVITATIONAL POLARONS 



A. Semiclassical description of fermions coupled to the gravitational field 

In the previous section, we considered the quantum system of Dirac electrons interacting via the attractive Coulomb 
potential. This system may appear in certain nonrelativistic models at the quantum phase transition between the phases 
of the fermionic models with different values of topological invariants 15511 . However, the Coulomb interaction breaks 
the emergent Lorentz symmetry. It is interesting, therefore, to consider the relativistic extension of the model defined 
by the partition function ( IIII.3I ). One such extension is discussed in this section. 

We consider the relativistic Dirac fermion interacting with the gravitational field. The action of a Dirac spinor in 
Riemann space has the form 

S f = J (i^D^)-mi>'il))\E\d 4l x. (IV.l) 

Here \E\ = det-E^, where E® is the inverse vierbein, 7^ = E£-f a . The covariant derivative is = 8^ + jwjf 7[ a 7b] , 
with 7[ a 7(,] = 5 (7 a 7b — Ibla)- The torsion-free spin connection is denoted by u^. It is related to E® and the affine 
connection as follows: 

v v e% = d v Ei - t^e; + ^ hv E\ = 0, 

D [v Efo = d [v E^ ] + U % [v E b ri =0. (IV.2) 

This results in: 

Uabv = ^{Cabc - C cab + C bca )E^. (IV.3) 

Here c aoc = r) a dE%E v c d\ y E'K, g^ = E a ^E h v r] ah , and r? v — V?^ = 0; indices are lowered and lifted with the aid of g 
and E. 

The partition function of the model is given by 

Z = J dtp dtp dE e* Id ' x (im^-m^-^RlEl) ; 

where = j^D^. 

Remark IV.l. The integral over Grassmann variables ijj in continuum field theory is a subject of the additional dis- 
cussion. There are several ways to define the functional integral: via lattice discretization, via reexpressing it as a 
functional determinant, etc. In all these cases, the presence of nontrivial metric leads to additional difficulties. Below 
we assume that the functional integral is defined in such a way that 



dtp d^ e< fdix l^+W = Det Q = U n X. 



(IV.4) 
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where A„ are eigenvalues of Q. (The spectrum of Q is discrete if we consider the system in the finite four-volume 
V4 = Jd 4 x \E\.) We imply the toroidal topology and the following boundary conditions in a synchronous reference 
frame: antisymmetric in time and symmetric in space coordinates. Eq. ( IIV.4b can be rewritten as 

Det Q = f # dip e* A » S d * x 1*1** (*)*»(*)e»c» 



dcdcDet^^W^" 5 "^, (IV.5) 
d(ip,ip) 

where we used the decompositions 

^ = }Z *n(x)Cn, *P + =Y, ^n(x)Cn. (IV.6) 

n n 

Here ^ n {x) are eigenfunctions of Q corresponding to eigenvalues A„, while c„,c„ are new Grassmann variables. 
Operator Q is implied to be hermitian with respect to inner product (t/j, (f>) = J d A x \E\ip + cj). Therefore, the eigen- 
functions satisfy J d A x \E\^^(x)^f m (x) = S nm . The key supposition about the integration measure over tjj is that 

with this normalization one has Det §^jk = 1- Eq. ( IIV.5b also allows us to calculate various correlation functions 

(ip + (xi) . . . i/j(xn))- Namely, we first represent ip as the series (IIV.6I 1. and then the integral over c, c is evaluated as 
usually. 

Again, the integral in the Euclidean space is not convergent, and should be redefined at high energies; see the 
discussion above. In order to bring the theory into the form suitable for the considerations similar to that of ll22ll . let 
us consider the system in the gauge corresponding to a synchronous reference frame. In this gauge, E®E?r] = 6°^; 
we also set E® = S® via the rotation of the reper in its internal space and the corresponding 5*0(3, 1) transformation 
of spinors. That is why the gauge is fixed both with respect to general coordinate transformations and with respect to 
the inner S0(3, 1) rotations of the reference frame. We denote 

7°[i^7°(a M + j<7[a7&]) ~m}= ida - H, 



-H = iEi^dj + i££ 7 V 2<7[c7 d ] - 7 °m, 



with the summation over j = 1, 2, 3. 

It is worth mentioning that the operator ido — His hermitian (while W is not). We have 



Z = JdE Det(i<9 -H)e-T^o *\ E \ 



In order to calculate the determinant Det(i<9o — T~L), we use anti-periodic in time boundary conditions. Suppose that 
we find the solution ( of the equation (ido — %)C = sucn mat Cu(t + T) = e~ inT C,n- (Here D.T is the Floquet index 
||22|l ). Then * fej o = e*^ (2k+x)t+ifit£ n [ s t h e eigenfunction of the operator (id - H): 

{ido - n)V k ,n = - (J(2A + 1) + fi)* fcl n. (IV.7) 

This results in 

Det(id Q -H) = n n cos(f2f T/2)n fe J(2fc + 1). (IV.8) 
Partition function takes the following form: 

Z = const fdEcxp(-im% fd i x\E\R + i^J2 n n - iT /2 Knn ™) e ~ e ° T (IV9) 

{K„}=0,1^ J n n 
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Here 80 is the constant that depends neither on the gravitational field nor on T. Following I22I1 . we interpret Eq. ( IIV.9b 
as follows. Numbers K n represent the number of occupied states with the Floquet index Q^T. These numbers may 
be 0, 1. The vacuum here corresponds to the negative "energies" (Floquet indices) occupied and positive "energies" 
empty. 

The physical meaning of the numbers K„ here is the same as in the previous section. The only difference is that 
now the Floquet indices appear in place of the energy levels and that the gravitational field depends on time. The 
semiclassical approximation for the gravitational field now leads to the variational problem 

= 5 {-Yy K n-\]^nT-m 2 P fd 4 xR\E\) 

n 

= 5 /^{E[ X ™-^*M2-( o- H + ^( 2A: + 1 ))*^ -m 2 P R}\E\. (IV.10) 

n 

Here E is varied, k is arbitrary, the normalization is J d A x \E\^~^ k q^e = T. Let us also introduce the wave 
function 

e -i(*(2*+i)+A B )t $Mjjf = CnAB . (IV.ll) 
The values A„ play the role of Lagrange multipliers. At Sl„ = A„, functions ( satisfy the following conditions: 

Jd i x\E\(E°^( = T, 
[0 - m}( = 0. 

Eq. ( II V. 1 Ob can be rewritten as 

= fd^fj^iKn-Wn^M^-H^Kn^E-mliSRlEiyY (IV.12) 

J n 

One can see that the variation of Q does not enter this expression. Therefore, for the determination of both Q^e and E 
we may use the variational problem 

= 6 [d 4 x{J2(Kn-hUi]P-™]Cn 
J n 

-m%R}\E\ (IV. 13) 

where the gravitational field and the wave functions £„ are varied. The wave functions are normalized so that 
/ d A x |i?|C + C = T. The additional constraint is that the wave functions are different, so that there are no states that 
are occupied more than once. The variation is performed with the fixed numbers K n . The final form of the variational 
problem is gauge invariant, although it was derived in a synchronous reference frame. As a result the gravitational 
field is defined by the Einstein equations with the energy-momentum tensor defined by the set of one-fermion states 
that represent the sea of occupied energy levels and the fermion-antifermion excitations given by the set K n . Fermion 
wave functions are defined by the Dirac equation in the given external gravitational field. 

After the variational problem is solved one says that the state with the found wave function ( n is occupied if K n 7^ 
for the corresponding value of n. In Eq. (|IV.91 > we need to sum up all such configurations with different arrays K n . 
The calculated values of E are to be substituted into the exponent in Eq. flIV.9b while the values of f2^ are given by 
Vt^ = C+'HCn- This is the generalization of the Hartree - Fock approximation. 

We came to the following result: 
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Theorem IV. 1. The partition function for the system of Dirac fermions interacting with the gravitational field is given 
by Eq. AIV.9\) . In the semiclassical approximation (for the energies E -C rap) the stationary phase approximation may 
be applied, leading to the variational problem Eq. (IIV.13\) . 

Remark IV.2. The variational problem Eq. dlV. 13b is gauge invariant, while the normalization of the wave function is 
not. We need 

Jd 4 x\E\C + C = T, (IV. 14) 

where the spinor £ and the time extent T are defined in a synchronous reference frame. Here T is a global characteristic 
of the space-time, |i?|c? 4 x is the invariant 4-volume, while C + C — C^al a C i s the time-component of the 4-vector. 



B. Gravitational one-polaron problem 



In order to investigate one-polaron state, we consider the two-point Green function 



G(t 2 -*i) = ^ J d$ d^P dE e l J^IBlW-mW-nto fd*x R\E\^+ ( tlyX )^(t 2 ,x) d 3 x\E(t 2 ,x)\. (IV.15) 

Here again the system is considered in a synchronous reference frame. It is implied that the corresponding terms are 
present in the integration measure over E. With all the notations introduced above, we can represent G as follows: 

Z q {K„}=0,lV 

where 

F E (h,t 2 ) = - [ [^(t u x)]+^(t 2 ,x)\E(t 2 ,x)\d 3 x 7 



Z = 



E 



dip e 



, f^x \E\R-iTY, n (K„-±)n* 



Here large value of mp allows to calculate integrals over E in the stationary phase approximation. When varying 
the terms in exponent, it is necessary to take into account F E (ti, t 2 ) which disturbs the effective action at t = ti,t 2 . 
However, if t 2 — ti = T, then, as in the previous section, we come to the following simplification: 



G ( T ) = \ E E / d <f ex p(- im P J d A x \E\R - iT^2[K n - 1/2]Q, E - iTD,^. (IV. 



16) 



q {K n }=0,l n^q 

The lemma follows: 



Lemma IV.3. In the quenched approximation, at the energies much less than m p , the one-polaron problem is reduced 
to the variational problem 

= 5 J "d i x{({ilp-m]C-mpR}\E\. (IV.17) 
Here C is the fermion wave function normalized according to Remark \lV.2\ 
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C. Newtonian limit 

In the non-relativistic limit, the energy-momentum tensor for the Dirac field has the form T^ v = C,i^ ,J 'd v 'C, 
mC^ Q8^ 8 vQ . The gravitational field is considered in the linear approximation g^ v = 77^ + f^, where r]^ 
diag (1, -1, -1, -1). Eq. ( IIV.17b in the gauge d^h^ = (where = - ^V^fp) has the form: 



= 8 J d*x {C[i0 - m - - ^[(d^U) 2 - \(d,f») 2 }}. (IV. 18) 

Here the variation is performed with respect to the wave functions Q and with respect to the graviton cloud /. As a 
result the graviton cloud is formed in accordance with the (linearized) Einstein equations with T^ v = (ij^d^^. 
In the non-relativistic limit (we neglect gravitational waves that are not caused by the given spinor field): 

/ o ° = 20, fZ = -2<f>8l a, b=l, 2, 3; 

iF u 7{^ } C ~ 20m 7 °C 
We arrive at the following system of equations: 

[i0-m- m4rf]C, = (IV. 19) 

Next, we require that <fr does not depend on time and denote m(f> = etp, e = As a result we arrive at 

/ idtC = -ia. ■ VC + mf3( + e<p(x, t)((x, t) m/om 

\A<p(x,t)=eC(x,t)((x,t) (1V - 2U) 

with the normalization 

f d 4 x CVCv^i = t. (IV.21) 

Here £' is spinor field in the synchronous reference frame, 7* = E®j a is the time component of covariant gamma- 
matrices (also, in the synchronous reference frame). 

The normalization of the spinor field is the subject of careful investigation. In the reference frame defined by 
harmonic gauge, both the spinor field and the gravitational field <fi are independent of time. However, in general case 



this is not the case in a synchronous reference frame. We may represent C'7*C = Cl^Q^oQr = J 1 * d gtd ■ wner e the 
current J M = C7 M C is defined in the original reference frame. In the weak coupling, we have in the original reference 

frame ^g ~ 1 - 2(j> and J^^P ~ C + C 

The latter follows from the Hamilton-Jacobi equation that defines the synchronous reference frame 

g^ v d ^ x l d j^ x l = 1. In the Newtonian approximation (f> « 1, and 9 ^ x } , v = 1,2,3 are of the same order as 
4>. Therefore, up to the terms linear in 4> we get 3°°[9 [a;'] ] 2 = [^o^o^'] ] 2 = l - A l so J' 1 if 1 = 1, 2, 3 are of the 
same order as <f>. Therefore, up to the terms linear in 4> we have J M = J° d j£j> = (^(E^dolx'] ~ £ + C 
That is why we come to the following normalization (valid in original reference frame in the weak coupling): 



1" 



/ d?xCCV~9 a [d 3 x( + ai - —<f) = 1. 
J J m P 



(IV.22) 



Lemma IV.4. In the Newtonian limit in harmonic gauge, the gravitational polaron problem is reduced to the system 
of equations MII.131 MIL 141 MV22\ with e = 
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V. STABILITY OF SOLITARY WAVES IN THE DIRAC-COULOMB SYSTEM 



A. Existence of solitary waves 



In this section we rescale fields in such a way that e = 1 (this is achieved when we substitute £ by ^( and ip by ^ip). 
As a result instead of the normalization condition Eq. ( MI. 15b (polarons in condensed matter systems) we have 



d 3 x( + ( = e 2 



For gravitational polarons we would have the constraint 



d 3 x( + ({l <p) =e z 



emp 



(V.l) 



(V.2) 



We express tp = A 1 \(\ 2 , where A 1 is the operator of convolution with — A ^ x ^ , and write the Dirac-Coulomb 
system dll-2b as the following Dirac-Choquard equation: 



idtC = -ia ■ VC + mf3( + ((x, ^A -1 ^, t)\ 2 , 



(V.3) 



with x £ M 3 , ((x, t) G C 4 . The solitary wave solutions 4> UJ e~ tu ' t with u> < m could be constructed by rescaling from 
the solutions to the nonrelativistic limit of the model. Such a method was employed in l29H37ll for the nonlinear Dirac 
equation and in l43l l44ll48tl for Einstein-Dirac and Einstein-Dirac-Maxwell systems. Let us mention that the solitary 
wave solutions to dV.3b with cj < m correspond to the solitary wave solutions of Dirac-Maxwell system with u> > — m 
when the magnetic field is neglected; such solitary waves were numerically obtained in j3~5ll . The sign change of lo 
is due to the different sign of the self-interaction (in the Dirac-Maxwell system, the self-interaction is repulsive for 
u> < m and attractive for u > — m; in the Dirac-Choquard equation dV.3l ), it is the opposite). The profile of the solitary 
wave C(x,t) = (f) U j(x)e~ luJt satisfies 



= -i oPdjfpu + m/3(t) u + <f>u,A 1 \(f> u 
j'=i 



(V.4) 



Let (j)^ (x) 



<f> p (x,uj) 



with (f> e , <p p £ C . In terms of <j> e and <f> p , JV.4t is written as 

w(j>e = -ia- ■ V</>p + m<f) e + faA' 1 (\<j) e \ 2 + \<t> p \ 2 ), 
uj(j) p = -icr ■ V4> e - m4> p + p A _1 (|0 e | 2 + \4> P { 2 ), 

with a ■ V = Y2j=i a fih w i m a o me P au li matrices. Let e > be such that e 2 = m 2 — oj 2 . We introduce functions 

# e (y, e), $ P {y, e) € C 2 by the relations 

(f> e {x,uj) = e 2 @ e (ex,e) 1 4> p {x,u) = e 3 <P p (ex, e); 

where e = V m 2 — w 2 . Let <9 yJ , A. y be the partial derivatives and the Laplacian with respect to the coordinates y = ex, 
so that d x j = ed y j, A x = e 2 A y . Then equations on <fr e , <fr p take the form 



m + to 



-ia ■ W y <Pp + # e A^ 1 (|#eI 2 + e 2 \<P p \ 2 ), 



(V.5) 



(m + u)$p = -icr ■ V v <P e + e 2 $pA-\\$ e \ 2 + e 2 \$ p \ 2 ). 



(V.6) 



Let u £ H 



00 fm3 



be a spherically symmetric strictly positive smooth solution to the Choquard-Pekar equation, 
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Such a solution exists due to JlJ, Qj, |34J] . Pick a unit vector n e C 2 . Then 

l e = TOe ir°(R 3 ,c 4 ), 



1 

2m' 



: 4 ) 



is a solution to dV.5l ). dV.6t corresponding to e = 0. The perturbation theory as in the above-mentioned references 
allows to construct solutions to (IV. 4b with ui € (cjq, to), with some cjo < to, such that 



4>e{x,Uj) 




(j) p {x,Uj)_ 





e 2 £ e (ex) + o(e 2 ) 
e 3 ^ p (ea;) + o(e 3 ) 



(V.8) 



where w and e are related by w = \Jm 2 — e 2 . 

Remark V. 1. Let us mention that among the solitary waves considered above, only those with the discrete values 
0J n ,K,e € (wo, to) satisfy constraint Eq. ( IV. j} (or constraint Eq. ( IV.2I )). These values are parametrized by the number 
of nodes n of the corresponding solution to the Choquard-Pekar equation ( 1V.7I ), the quantum number k = ±1, and 
the value of charge e that enters Eq. dV.lt (or constraint Eq. dV.2b ). This pattern is described in detail in ll3~5ll (in the 
context of the Dirac-Maxwell system). 



B. Linear stability of solitary waves 

We assume that cjo < to is such that for uj E (uio, m) there are solitary wave solutions to dV.3t . <p cu (x)e~' lu ' t . Taking 
the Ansatz ((x, t) = (cf) u (x) + p(x, t))e~ lult we derive the linearization at the solitary wave <f) UJ (x)e^ lult : 

ip = (D m -cj + A-^l^l 2 ))^ + A^Go*^ + 

where D m = —ia ■ V + to/?. We are looking for the eigenvalues of the linearization operator in the right-hand side. 
That is, we substitute p(x, t) = £,(x)e xt , with £ € L 2 (M 3 , C 4 ), £^ 0, getting 

iA£= (D m - w + A- 1 (|0 w | 2 ))e + ^A- 1 (r^ + CO, (V.9) 

and we would like to know possible values of A. If there is Re A > corresponding to £ ^ 0, then the linearization 
at a solitary wave is linearly unstable, and we would expect that the solitary wave is ("dynamically") unstable under 
perturbations of the initial data. 

Theorem V.l. There is u>\ £ [uio, to) such that the "no-node" solitary waves with lo £ (cJi, m) are spectrally stable, 
so that in dV.9t one has Re A = ( unless £(x) = 0). 

Remark V.2. By the "no-node" solitary waves we mean the solutions dV8t constructed from the strictly positive 
solution to the Choquard-Pekar equation dV.7t . 

Let us mention that such "no-node" solutions that satisfy constraint Eq. d V. 1 b (or constraint Eq. ( IV.2b ) only exist if 
the value of e 2 is sufficiently small. That is why the above results on existence and stability may be reformulated as 
follows: 

Theorem V.2. There is > such that the Dirac-Coulomb system dll.2t with e 2 £ (0, eg) has "no-node" solitary 
wave solutions £(x, t) — <j){x)e~ lu}t which satisfy the constraint (or AV.2j ) and are spectrally stable. 

Remark V.3. According to the scaling in the Ansatz dV.8t , one has J \<j> u | 2 d 3 x ~ e ~ (to — uj) 1 / 2 , hence e 2 and w in 
Theorem lV.2l are related by 

e 2 ~ (to — uj) 1 / 2 , uj < to. 

We point out that dV9t is R-linear but not C-linear, because of the presence of £* . Let us rewrite dV.91 ) in the C-linear 
form. For this, we introduce the following notations: 

h 
-h ' 



1 = 



ReC 



Re< 
Im < 
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Re aP —Tm a 3 
1hq j Re a 3 



Re/3 -Jm0 
Jm0 Re/3 



Then dV.91 > can be written as 



where 



XI, = JL(w)t, 



(V.10) 



L(w)t = (D m -w + A~ 1 (\$J 2 ))t + 2A- 1 ((b^)4)., 



3 

D m = ^ Jo 3 ' 3, + pm. 

3=1 

The operators D m and L(cj) are self-adjoint. 

Theorem lV.ll is the immediate consequence of the following lemma. 

Lemma V.4. Let Wfe G (0, to), fc G N; Wfc — > to. There is no sequence Afe £ er p (jL(u;fc)) vv/f/z Re Afe 7^ 0. 

Proof. The lemma is proved in several steps, which we now sketch; more details will appear in ifTHl . 

First, one shows that if there is a sequence of eigenvalues A& G a p ( JL(wfe)) such that lim Afe exists, then 

k— >oo 

lim A fe C {0,±2mi}. 

k^oo 

The proof of this statement follows from the fact that in the limit u — > to, as ||^>a;||z,oo — > (Cf. (IV.8I 1). the operator 
jL(w) turns into J(D m — to). According to [53], there is the limiting absorption principle for the free Dirac operator 
D m : its resolvent, (D m — z)~ l , is uniformly bounded in certain weighted spaces, uniformly for | Kbz\ > to + S (for 
any fixed 6 > 0) and Im z 7^ 0. This implies that the resolvent of JL(w) is bounded in these weighted spaces outside 
of the union of iM. with open neighborhoods of "thresholds" A = and A = ±2mi, as long as ui is sufficiently close 
to to. In turn, this implies that as ojk — > to, the eigenvalues Afe can not accumulate but to these three threshold points. 

Further, the eigenvalues Afe with Re Afe 7^ can not accumulate to ±2mi. This follows from the fact that if 
Afc — > Ab € iM\0 as Wfc — >• G iK\0, then Ab itself has to belong to the point spectrum of JL(w(,) (corresponds to 
the 1? eigenfunction); this result is again based on the limiting absorption principle. At the same time, there can be 
no L 2 eigenfunctions of a constant coefficient operator J(D m — to). 

Finally, one has to study the most involved case Afe — > 0. One first proves that if Afe — > and Re Afe 7^ as uik — >• to, 
then necessarily Afe = 0(m — ujk)- Then one studies the rescaled equation. The conclusion is that the families of 
eigenvalues for the linearization of the Dirac-Choquard equation, Afe G <7 p (JL(cjfe)), are deformations of families of 
eigenvalues for the linearization of the Choquard-Pekar equation (which is a nonrelativistic limit of Dirac-Choquard 
equation). The presence of eigenvalues with nonzero real part in the linearization of Choquard-Pekar is controlled 
by the Vakhitov-Kolokolov stability criterion lf50ll : for the linearization at no-node solutions, this criterion prohibits 
existence of such eigenvalues. This finishes the proof of the lemma. □ 

We reproduce the Vakhitov-Kolokolov stability criterion f50ll in application to the Choquard-Pekar equation in the 
Appendix (see Lemma lATI below). 



VI. CONCLUSIONS 



In the present paper we considered solitary waves in the system of Dirac fermions interacting via the Coulomb 
attraction from two points of view: the point of view of a physicist and the point of view of a mathematician. 

On the physical side the solitary waves describe polarons that may appear in two situations. The first one corre- 
sponds to certain condensed matter systems when massive Dirac fermions emerge that interact with optical phonons. 
Also the polarons in the system of true relativistic Dirac fermions interacting with gravity in the Newtonian limit are 
described by the above-mentioned solitary waves. 

On the mathematical side we develop analytical methods for the investigation of solitary waves. These methods 
are based on the observation that these soliton solutions are obtained as a bifurcation from the solitary waves of the 
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Choquard-Pekar equation. Basing on this approach, we demonstrate that the no-node gap solitons for sufficiently 
small values of e are spectrally stable. 

It is worth mentioning that the solitary waves similar to the considered in the present paper may also exist in two- 
dimensional systems like the boundary of the topological insulators or the graphene. This may occur if the interaction 
between the electrons of the 2D system with the balk phonons (with the substrate phonons in the case of graphene) is 
strong enough. We postpone the consideration of the corresponding 2D solitary waves to future publications. 

The authors would like to kindly acknowledge the private communication with G.E. Volovik, who prompted to 
consider the relation of the polaron problem to the solitary waves in Dirac-Coulomb system. A.C. greatly benefitted 
from discussions with Gregory Berkolaiko, Nabile Boussai'd, and Boris Vainberg. The work of M.A.Z. was partially 
supported by RFBR grant 1 1-02-01227, by the Federal Special-Purpose Programme "Human Capital" of the Russian 
Ministry of Science and Education, and by the Federal Special-Purpose Programme 07.514.12.4028. 



Appendix A: Vakhitov-Kolokolov criterion for the Choquard equation 

We consider the Choquard-Pekar equation, 

^ f C = -^AC + mC + CA- 1 |C| 2 , 
2m 

where C( x it) <= C 3 and x <G M 3 . We are interested in the solitary wave solutions CO;*) = u ul {x)e~ 
satisfies 

- (m - uj)u u = -—Auu + A _1 (m^ )u u . 
Im 

Given a solution to ( I V.7b . then, for any uj < m, the profiles 



(A.l) 



uj e K; u u 



(A.2) 



u u (x) = (m — u))y/2mu(\x\\/2m(m — uj) ) (A. 3) 

correspond to a family of solitary wave solutions to iA.2i , Note that this scaling is the same as that of <f> e in ( IV.8I ). 



Lemma A.l. For uj < m, the no-node solitary wave solutions C, u {x, t) = u u (x)e lu)t to the Choquard-Pekar system 
are spectrally stable. 

The linear stability of no-node solitary waves of the Choquard-Pekar equation follows from the Vakhitov-Kolokolov 
stability criterion IdOTi which is applicable to systems of the Schrodinger type. It also follows from fPjll (where the 
orbital stability of these solitary waves is proved). Let us sketch the argument. First, we notice that, by ( |A.3t , the 
charge of the solitary wave u UJ (x)e~ lult is given by 







(x)| 2 d^x ~ (m — uj) 1 / 2 , uj < m; 



therefore, 



dQ{u M ) 
duj 



< 0, 



uj < m. 



(A.4) 



We consider the solution to the Choquard-Pekar equation in the form of a perturbed solitary wave, 

CO, t) = {u u {x) + R(x, t) + iS(x, t))e- w \ 
with R, S real-valued. The linearized equation on R, S is given by 



JIM 



where 



1 = 



1 

-1 



Li(w) 
Lo(w) 



(A.5) 



(A.6) 
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and 



L (u) 



2m 



+ A- 1 (^), 



Li(w) = L (ui) + 2A l {u u -)u u . 

Both operators Lq(uj) and Li(oj) are self-adjoint, with <r ess (Lo(u)) = <r ess (Li(oj)) = [m — u,+oo). Clearly, 
Lq(w)u u1 = 0, with € Od(L ) an eigenvalue corresponding to a positive eigenfunction u u ; it follows that is a 
simple eigenvalue of L , with the rest of the spectrum separated from zero. Taking the derivatives of the equality 
Lq(w)uu = with respect to Xj and uj, we get: 

L!{uj)d Xj u w = 0, Li((jj)d u u u = u u . (A.7) 

The first relation shows that Ai = is the point eigenvalue of L\(u), and since d Xj u u vanishes on a hyperplane 
Xj = 0, there is one negative eigenvalue Ao < of Li(u>), 



Now we may determine the spectrum of j\(u) = 



L (w) 




We closely follow |50D. If 



is an eigen- 



function corresponding to the eigenvalue A G C, then —\ R = LqL\R. If A ^ 0, then one concludes that R is 
orthogonal to kcr Lq = Span(ii u ), hence we can apply Lq ; taking then the inner product with u u , we have: 



X 2 (R,Lq 1 R) = (R,LxR). (A.8) 



With Lq, Li being self-adjoint, this relation implies that A € R. Since Lq(uj) is non-negative and R _L kerLo> one 
has (iiu, L^uJ) > 0. The solution to 

H :=mt{{R,Li{L0),R); \\R\\ = 1, =0} 

satisfies L\{ui)R = //i? + vu^, where \i, v S R play the role of the Lagrange multipliers. Due to the condition 
(u u , R) = 0, /i delivers the zero value to the function 

/(z) = (u u , (Li(u) - zy 1 ^}, z e 

with p(Li) denoting the resolvent set of L\. Since kcrii is spanned by djU^ and therefore is orthogonal to u u , we 
can extend f(z) to z G (Ao, A2), where Ao = inf a(Li(ui)) < and A2 is the first positive eigenvalue of L\ (or the 
edge of the essential spectrum, A = m — uS). We need to know whether fi is positive or negative. Since f'(z) > 0, the 
sign of fj, is opposite to 

/(0) = = {u u ,d u u u ) = ^^j"" . 



In the second equality, we used the second relation from (IA.7I ). From (|A.4t , we conclude that /(0) < 0. Thus, fi > 0. 
By (EOl A 2 < 0, leading to C iR. 

This shows that there are no families of eigenvalues of with nonzero real part bifurcating from A = at 
uj = in. Since bifurcations of eigenvalues from A = for the linearizations of the Choquard-Pekar equation and 
the Dirac-Choquard equation JV.31 > (which is equivalent to the Dirac-Coulomb system) have the same asymptotics as 
uj — > m, we conclude that neither are there families of eigenvalues of JL(u;) with nonzero real part. 

Remark A.2. The rigorous proof of linear stability of solitary wave solutions to the Dirac-Choquard equation requires a 
more detailed analysis of the spectrum of Namely, one needs to know whether there are resonances or embedded 
eigenvalues. Theoretically, resonances or embedded eigenvalues of higher algebraic multiplicity could bifurcate off 
the imaginary axis into the complex domain, yielding a family of eigenvalues Afe € er p (JL(o;fc)) with Uk — > m, 
Afe = 0(m — cuk), Re Afc 7^ 0, (and resulting in the instability of Dirac-Choquard solitary waves), although we expect 
that generically this does not happen. 
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